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INTRODUCTION 


A computer  simulation  of  a two  rate  gyro  stabilized  platform,  laser 
semi-active  seeker  Is  given  here.  The  seeker,  shown  In  Figures  1 and  2, 
Is  the  Wide  Fie!d-of-V!ew  (WFOV)  seeker  developed  for  use  In  the  Terminal 
Homing  Accuracy  Demonstration  (THAD)  program.  A thirteenth  order  seeker 
dynamic  model  Is  studied  and  compared  to  actual  seeker  performance.  The 

method  of  Chen  and  Shleh^  Is  used  to  simplify  this  model  to  a second 
order  system.  The  simplified  model  was  then  compared  to:  (1)  the 

thirteenth  order  model,  (2)  a contractor  third  order  model,  and  (3) 
the  actual  system  response.  The  agreement  between  the  responses  was 
quite  good  and  Is  discussed. 


This  report  emphasizes  the  digital  programs  used  In  this  study 
since  they  can  be  used  for  a wide  variety  of  engineering  applications. 
These  routines  Include  a Runge-Kutta  numerical  analysis  method,  a 
plotter  package,  a compound  polynomial  simplification  process  and  a 
transfer  function  order  reduction  program.  A brief  description  of  each 
program  follows  with  a sample  problem  Illustration  where  appropriate. 

2.  RUNGE-KUTTA  DIFFERENTIAL  EQUATION  SOLUTION 


To  obtain  the  seeker  response  In  the  track  mode,  a fourth  order 
Runge-Kutta  algorithm  for  the  thirteenth-order  servo  system  was  used. 
Figure  3 shows  the  differential  equation  solution  program  for  the 
seeker.  All  Initial  conditions  of  the  state  variables  XN(1)  through 
XN (13)  were  set  to  zero  radians.  The  program  was  set  up  to  print  values 
of  the  state  variable  In  degrees  every  10  msec  using  a Runge-Kutta  time 
-A 

Increment  of  h « 10  second. 

To  measure  the  seeker's  tracking  rate,  a laser  radiation  source  Is 
set  In  a position  which  Is  3°  from  the  seeker's  caged  1 tne-of-slght.  In 
the  simulation  this  angular  error  (A10S  * 3®)  Is  Initially  read  Into  the 
program.  Then  at  time  t ■ 0 the  seeter  Is  switched  Into  the  track  mode 
permitting  a track  command  signal  (SAMP  ■ £ 2v)  to  align  the  sensor  head 
with  the  source.  At  t « A5  ms,  SAMP  « 0 for  a period  of  5 ms  In  a 
manner  Identical  to  the  seeker's  signal  processor.  Then  at  t “ 50  ms 
the  seeker  compares  the  ALOS  error  to  the  platform  position  XH 1(1)  to 
determine  whether  SAMP  will  be  positive  or  negative.  This  comparison 
continues  to  occur  at  50  ms  Intervals  to  update  the  seeker  head  posi- 
tion. 


{ 1 ] C.  F.  Chen  and  l.  S,  Shleh,  "A  Novel  Approach  to  Linear  Model 
Simplification,"  Int.  J.  Control,  Vol.  8,  1968,  pp.  561-570. 
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The  state  equations  FX ( 1 ) through  FX ( T 3)  are  taken  from  Figure  2, 
and  Initially  set  equal  to  zero.  The  state  equations  are  defined  as: 

XO)  » x 

FX ( 1 ) = dx/dt 


FX (12)  = d12x/dt12 
FX  (13)  - d13x/dt13 


Drift  rates  are  entered  simply  by  adding  a term  to  the  twelfth 
state  equation  such  that 

FX (12)  = 1.E5/3.77  * (-8.58E-3  * x(12)  X(11)  + x(2)  + DRIFT) 


Similarly  by  making  changes  In  the  state  equations  In  accordance  with 
the  constants  of  Figure  1,  either  the  pitch  or  yaw  response  can  be 
observed.  The  pitch  response  will  be  discussed  throughout  this  report 
although  both  cases  were  studied  and  found  to  be  essentially  Identical. 

Data  cards  are  read  Into  the  programs  as  follows: 

1.  IPPS  - Pulses  per  second  rate  of  laser  Illuminator  (IPPS  » 

20). 

2.  N - Order  of  differential  equation  (N  • 13). 

3.  AtOS  - Angular  1 Ine-of-slght  error  In  degrees  (ALOS  * 3.0). 

4.  TN,  H,  XN(I)  ...  XN 03)  - Initial  time  (TN  - 0),  Runge-Kutta 

Increment  (H  * 10  **),  In  1 1 la ' state  variable  conditions  In 
radians  (XN(I)  ...  XN (13)  are  all  0). 

These  values  are  printed  by  the  program  for  reference,  as  well  as 
the  calculated  sampling  period  (TSAW  «*  .050  second)  for  SAHP  sign 
changes . 
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3. 


PLOT  PACKAGE 


An  easily  used  plot  package  Is  Included  In  the  Runge-Kutta  analysis 
program  to  provide  a graphic  recording  of  the  system  response  (Figure  4). 
This  routine  provides  a neatly  labeled  plot  of  up  to  five  functions  of 
501  points  each  versus  time,  complete  with  grid  lines  and  polnt-by-polnt 
print  out  of  the  first  two  functions  alongside  the  graph. 

To  use  the  plotter  simply  call  PLOT  (Y,  M,  NF,  NS). 

Y ■ the  array  of  data  to  be  plotted.  The  first  function's 
data  Is  In  the  form  Y(l,l),  Y (2 , 1 ) ...  Y(NF,1),  with  the 
second  function  (If  any  Is  to  be  plotted  on  the  same  graph) 
similarly  arranged  as  Y(l,2),  Y(2,2)  ...  Y(NF,2),  and  so 
on  with  as  many  as  five  functions. 

M » number  of  functions  per  graph  (up  to  5). 

NF  ■ number  of  points  per  function  beginning  with  the 
Initial  conditions  (up  to  501  po*"^). 

NS  ■ maximum  upper  limit  of  the  ordinate  axis. 

It  Is  Important  to  remember  that  this  routine  plots  from  NS  - 100  to 
NS  (e.g..  If  NS  ■ 50  the  ordinate  Is  labeled  from  -50  to  +50;  NS  * 100 
the  axis  Is  from  0 to  100). 

The  Initial  conditions  of  the  platform  output  XNl(l)  and  SAMP 
are  loaded  Into  the  data  array  elements  Y(!,l)  and  Y(l,2)  early  In  the 
Runge-Kutta  program  at  lines  53  and  54,  while  the  remaining  values  of 
these  two  functions  are  called  Into  the  plot  program  by  C(KK+1,I)  and 
C(KK+1,2)  at  lines  126  and  127. 

The  plot  parameters  used  In  the  Runge-Kutta  program  are  given  as, 

H » 2 for  the  two  functions  to  be  plotted  versus  time  on 
the  same  graph,  XNl(l)  and  SAHP  respectively. 

NF  « 401  for  a total  of  40!  points  per  function. 

NS  * 50  normally  labels  the  graph  from  -50  to  +50.  A 
discussion  of  the  use  of  NS  with  scaled  functions  will 
be  given  in  the  next  section. 

Some  computer  systems,  such  as  the  COC  6600,  may  not  properly 
execute  the  plot  program  unless  the  Y data  array  dimension  statements 
are  the  same  In  both  the  main  program  and  the  plot  subroutine.  Further- 
more, the  Y array  may  have  to  be  dimensioned  Y(5,50!)  rather  than 
Y (501 ,5)  depending  again  on  the  system  being  used.  Of  course  any 
desired  Increase  In  the  maximum  number  of  functions  or  points  to  be 
plotted  can  be  made  by  simply  Increasing  the  data  array  dimension 
statements  In  the  two  programs. 


i*.  USE  OF  THE  PLOT  ROUTINE  WITH  SCALE  FACTORS 


Any  functfon(s)  may  be  scaled  for  clearer  graphical  display  by  the 
plot  subroutine.  This  requires  the  simple  change  of  only  the  plot  data 
array  cards  and  an  ordinate  axis  labeling  card  within  the  subroutine. 

A description  of  the  steps  used  to  scale  the  platform  position  XN1(1) 
follows  which  serves  as  a guide  for  any  scaling  that  may  be  desired. 

(1)  Referring  to  lines  53  and  126  of  the  Runge-Kutta  program, 
the  platform  angular  position  was  multiplied  by  a scale  factor  of  10 
to  show  Its  relationship  to  SAMP  more  adequately.  Thus  the  data  array 
cards  for  the  Initial  condition  and  all  subsequent  values  of  XNl(l) 
become  respectively 


C(1,l)  - 10.  * XN 1 ( 1 ) 


and 


C(KK  + 1,1)  » 10.  * XN1  (1) 

(2)  Line  18  In  the  plot  routine  was  changed  to  label  the  ordinate 
axis  according  to  the  scaled  platform  function.  Normally  this  line  would 
read 


101  L(l)  - 10  * I - 110  + NS 

labeling  the  ordinate  axis  from  NS  - 100  to  NS.  However,  since  XN1(1) 
was  Increased  by  a factor  of  10,  the  axis  scale  had  to  be  reduced  by 
the  same  factor.  Thus  line  18  becomes 

101  L(l)  - (10  * 1 - 110  + NS)/10. 

which  now  labels  the  ordinate  from  (NS  - 100)/l0  to  N5/10.  Note  that 
NS  « 50  In  the  main  program  since  the  scaled  values  of  XNl(l)  range 
from  0°  to  approximately  30®.  The  change  In  line  18  of  the  plotter 
now  expands  the  ordinate  from  -5  to  +5. 

(3)  Line  52  In  the  plot  routine  rescales  XNl(l)  to  Its  original 
values  for  printing  along  side  the  graph.  Thus, 

Y(N,L)  - Y(N,!)/!0. 

Note  that  SAMP  was  not  scaled  to  a larger  value  and  therefore 
does  not  correspond  to  the  ordinate  label.  However  this  Is  unimportant 
since  only  the  sign  of  SAMP  rather  than  Its  magnitude  Is  of  Interest 
here. 
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5.  POLYNOMIAL  SIMPLIFICATION  ROUTINE 


It  Is  often  a long  and  tedious  process  to  obtain  the  total  poly- 
nomial transfer  function  of  omplex  systems.  Many  times  the  simpli- 
fication of  such  systems  Is  >**t  In  a form  which  is  a long  combination 
of  polynomial  sums  and  products.  The  following  program  (Figure  5)  Is  a 
useful  means  of  executing  the  addition  and  multiplication  of  these 
polynomials  In  order  to  produce  a transfer  function  of  the  form 


F(s) 


v"  + Vi5"'1  + + V + l!o 


Let  the  system  transfer  function  be  In  the  general  form  of  a quo- 
tient of  two  compound  polynomials 

(Aj  + Aj  + ...  + Am)  (B) 

F(s)  " Tc,  + 't2  cNJ  (6) 

where  A] , A2  ...  A^,  B,  Ci,  C2  ...  Cfj,  0 are  themselves  products  of 
polynomials  and  constant  factors.  For  example,  Aj  might  be  a product  of 
two  polynomials  and  two  constant  factors  such  that 


\ 2 

A1  * KlVal3S  + a12s  + aus  + al0^a2ls  + a20^ 

The  operation  of  this  program  Is  rather  straightforward.  First  the 
A polynomials  are  solved  (if  any  aro  present),  then  added  together  and 
finally  their  sum  Is  multiplied  times  the  B polynomial  product.  The 
resultant  numerator  Is  printed  giving  the  order  of  "s"  and  Its  coeffi- 
cient. The  denominator  polynomials  are  handled  In  the  same  manner. 

The  following  sample  problem  Illustrates  the  use  of  data  cards  for 
this  program. 

Let 

5(s*  - 2s3  + 2s2  + s - l)(s2  + 2) 

r(s)  » — -**->— — 1 « 

l(s4  + 2s  - 3) (s  + 2)  ♦ 20{sJ  - 2)(s  + l))(sz  - 3s  + 1) 
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where 


B “ 5(sA  - 2s3  + 2s2  - l)(s2  + 2) 

C,  - (s2  + 2s  - 3)(s  + 2) 

C2  - 20(s3  - 2) (s  + 1) 

0 ■ s“  - 3s  + l 

Data  cards  are  read  In  the  following  order  beginning  with  the  A polynomials 
If  any  are  present  (Figure  6). 


1. 

NA 

Number  of  A polynomials  which  must  be  added. 
Since  there  are  no  A polynomials  here,  NA  ■ 0 
and  we  now  consider  the  B polynomial. 

2. 

NP,  KK 

Number  of  polynomials  to  be  multiplied  In  B and 
the  number  of  constant  factors  respectively. 
Here,  NP  - 2 and  KK  « 1. 

3. 

K(l) 

Constant  factors  In  B (other  than  1 .0) . If  no 
constants  are  present,  simply  omit  this  data  card. 

4. 

Nl 

Order  of  the  first  polynomial  to  be  multiplied  In 
B.  Since  the  first  polynomial  Is  fourth  order, 

Nl  ■ 4. 

5. 

POLY  1(1) 

Coefficients  of  the  first  polynomial  factor  begin- 
ning with  the  highest  order  "s"  coefficient.  This 
card  reads:  1.E0  -2.E0  2.E0  1.E0  -1.E0). 

6. 

N2 

Order  of  second  polynomial  factor  In  B.  N2  - 2. 

7. 

POLY  2(1) 

Coefficients  of  the  second  polynomial.  This  card 
Is  1.E0  0.E0  2.E0. 

8. 

NA 

Number  of  C polynomials  In  the  denominator.  NA  - 2 

9. 

NP,  KK 

Refers  to  Cj.  Here  NP  » 2 and  KK  • 0. 

K(l) 

Since  KK  • 0,  this  card  Is  omitted. 

10. 

Nl 

Order  of  first  polynomial  to  be  multiplied  In  C,. 

Ni  - 2.  * 

n. 

POLY  1(1) 

Refers  to  Cj. 
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12. 

NZ 

Refers  to  Cj. 

13. 

POLY  2(1) 

Refers  to  Cj . 

14. 

NP,  KK 

Refers  to  Cj. 

15- 

19. 

Refers  to  C2- 

20. 

NP,  KK 

Here  NP  » 2 (Includes  the  C polynomial  sum  and 
0 factors).  KK  = 0. 

the 

K(l) 

Since  KK  » 0,  this  card  la  omitted. 

21. 

N2 

Order  of  the  first  D polynomial  factor.  POLY  1(1) 
Is  the  sum  of  the  C polynomials  which  Is  stored  In 
the  program.  No  cards  are  needed  for  POLY  l(l) 
following  the  addition  of  A or  C polynomials.  Here 
the  first  D polynomial  Is  second  order,  so  N2  = 2. 

22. 

POLY  2(1) 

Coefficients  of  the  first  D polynomial  factors. 
This  card  ts  1.E0  -3.E0  1.E0. 

Referring  to  the  seeker  block  diagram  It  Is  apparent  that  the  trans- 
fer function  up  to  the  platform  rate  output  Is  of  the  twelfth  order.  In 
some  cases  It  may  be  of  Interest  to  monitor  the  tracking  rate  of  the 
platform  as  well  as  Its  actual  position.  For  this  reason  a twelfth  order 
system  polynomial  was  obtained  using  this  program  and  the  data  cards  In 
Figure  7.  Finally,  It  became  desirable  to  reduce  this  rate  portion  of 
the  system  to  a second  order  approximation  model  for  easier  design  mani- 
pulation. The  following  program  was  used  to  this  effect. 

6.  REDUCTION  IN  OROER  FOR  A DIFFERENTIAL  EQUATION 

» L 

This  program  reduces  the  order  of  an  n order  differential 
aquation  by  expanding  Its  s-plane  polynomial  Into  a continued  fraction 
and  truncating  certain  of  Its  quotlentslU.  The  remaining  expansion 
quotients  are  then  used  to  write  all  the  lower  order  polynomials  beginning 
with  order  n - 

The  high  order  polynomial  Is  Initially  arranged  In  ascending  order, 
beginning  with  the  constant  terms. 

pn  + V + p?sZ  + •••  ♦ P sm 
F(s>  - ! * 5L. 

Qq  + QjS  + QgS^  ♦ ♦ Qnsn 

From  this  polynomial  the  program  calculates  and  prints  the  2n  continued 
fraction  terms  (e.g.,  a 12th  order  equation  has  24  terms). 


8 


Am, 


F (s) 
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1 
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H 
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1 


1 
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To  derive  an  approximation  of  order  tu,  only  the  first  2w  terms  are 
retained  for  a reversal  of  the  expansion  process  which  yields  the  lower 
order  polynomial.  For  a second  order  model  of  the  twelfth  order  seeker 
platform  system,  the  first  four  continued  fraction  terms  were  saved 
so  the: 


F(s) 


1 


34.474  + ! 

-6.8275  + ! 

s 

-3.057  + ~ 

20.732 

s 


i 


i 

>'• 

% 


i 

P 


’ V*. 

/ 

'’'•./‘kKi.'yTfr''' 


which  produced 


F(s) 


13.90As  + 432.71 
s2  + 1*15. 96s  + 14917.6 


The  program  Is  designed  to  print  the  terms  of  the  original  system 
and  all  of  Its  lower  order  approximations  In  ascending  order. 

Data  Is  eastly  read  Into  this  program  as  Figure  9 demonstrates. 

The  only  restriction  Is  that  the  denominator  Is  assumed  to  be  one  order 
higher  than  the  numerator  for  execution  of  the  continued  fraction  process. 
This  condition  Is  satisfied  by  the  Introduction  of  "dummy"  data  terms. 

For  Instance  the  WFOV  seeker  platform  polynomial  has  a 9*^  order  numerator 
and  a 12th  order  denominator  such  that 
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F<s) 


♦ * » 


po  + pis  + 


Qq  + Q^s  + . . . 


+ PgS 


+ ai2S 
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To  make  the  numerator  one  order  less  than  the  denominator  simply  enter 
zero  for  the  numerator  s^O  and  coefficients  so  that 


F (s) 


V 


Qq  + Q,s 


+ ...  + Pgs9  + Os10  + Os11 

nr 


T 


IT 


' v ^10s  + ^10S  + ^US  + ^12S 


T2 


Referring  to  Figure  9>  data  Is  entered  for  the  seeker  platform  as 
fol lows. 


1.  NN  - number  of  terms  In  the  numerator.  NN  «=  12. 

2.  A(l)  ...  A(NN)  - coefficients  of  the  numerator  In  ascending 
order  beginning  with  the  constant  term.  Note  that  A ( 1 1 ) 
and  A ( 1 2)  are  dummy  (zero)  terms. 

3.  B(l)  ...  B(ND)  - coefficients  of  the  denominator  In  ascending 
order  beginning  with  the  constant  term  ND  » NN  + 1 (reference 
line  19  of  the  Figure  8 program  listing). 

7.  RESULTS 


The  thirteenth  order  model  of  the  WFOV  seeker  yielded  an  accurate 
simulation  of  the  seeker's  time  response  to  a 3°  1 Ine-of-slght  error. 

As  expected,  the  model  demonstrated  a 3Vsec  tracking  rate  with  an  Input 
forcing  function  (SAMP)  correction  every  50  ms  (Figure  10).  The  state 
equation  set  up  for  the  seeker  Is  given  In  Figure  2 as  equations  FX ( 1 ) 

through  FX(13),  where  FX(l)  **  , FX(2)  - ^4  , etc. 

oc  dr 

A polynomial  transfer  function  was  then  derived  from  the  servo 
block  diagram  and  Is  given  as 


r(») 


«5  ♦ p,.»*  ♦ n.t5  ♦ p.»*  ♦ p.t’  * 


V 4 v * p?‘  * v * y * v * y * Pj»  ♦ p,»  ♦ p0  , 

""■  IS"1  1 (|C«M|‘"“»‘-i||p"i|“‘LM»W|,!'i»l|j*»r,|,«4iMWW|Ii..'«,w.«*»,yiiw|  .■yiiiMM.iw  .■  Mini  tmiiin  , m 

ii*  * in*  4 «io*  ♦ y9  ♦ y8  ♦ y7  ♦ y6  ♦ y5  ♦ y*  ♦ y3  ♦ <yl  ♦ y ♦ ^ * 
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where 


q12 

- 1 . 265 ( 1 O-23 

qll 

- 3.327(10'19 

q10 

- 2.084(10'15: 

P9  - 1.556(10"15) 

q9 

- 5.7O3O0'12: 

P8  - 3.340(10-11) 

q8 

- 7.968(1 o'9) 

Py  - 8.450(10-8) 

q7 

- 5.966 (10~8) 

P6  - 7.229 (10'5) 

% 

- 2.440(10'5) 

P5  - 2.424 (10~2) 

q5 

- 7.930(10~1) 

P4  - 4.811 

q4 

■ 1 .625 ( 1 02) 

P3  - 5.802(102) 

q3 

- 1.965(104) 

P2  - 4.495 dO4) 

q2 

- 1.493(106) 

P,  - 1.909d06) 

ql 

- 6.1 16 (107) 

P0  “ 3.178(107) 

q0 

- 1 .095  do9) 

The  twelfth  order  polynomial  represents  the  system  up  to  the  platform 

rate  output  while  the  final  i Integrator  provides  the  platform  position 
(reference  Figures  1 and  2). 


The  twelfth  order  polynomial  was  reduced  to  the  second  order  approx I- 
mat  ion 


F(s)  - 1,3.904s  ♦ *>32.71 

s2  + 415.96s  + 14917.6 

and  Is  shown  In  Figure  11  as  part  of  a third  order  system  model.  The 
state  equations  for  Runge-Kutta  Integration  of  this  system  are 


II 


FX(3)  = SAMP  - 415.96  X(3)  - 14917.6  X(2) 

FX(2)  - X(3) 

FX(l)  = 13.904  X(3)  + 432.71  X (2) 

The  third  order  approximation  accurately  demonstrated  a 3°/sec 
tracking  rate,  however  Its  forced  platform  sweep  generally  remained 
below  a 3°  position  and  was  equal  In  magnitude  but  180°  out  of  phase 
from  the  higher  order  system  (Figure  10). 

The  discrepancy  between  the  two  models  was  due  to  the  slight  lead 
of  the  third  order  system.  For  Instance  at  time  t = 1.1  second  the 
platform  positions  of  the  thirteenth  and  third  order  systems  wern 
2.99992  and  3.00528  respectively.  Consequently  SAMP  remained  positive 
for  the  first  system  but  changed  to  a negative  value  for  the  latter 
system  causing  it  to  decrease  in  value.  This  out  of  phase  relationship 
maintained  itself  throughout  the  entire  forced  switching  mode. 

Oispite  the  above-mentioned  minor  differences  the  lower  order  model 
as  a whole  proved  to  be  a sufficiently  accurate  representation  of  the 
thirteenth  order  system.  This  confirms  the  validity  of  the  continued 
fraction  expansion  process  used  in  the  fourth  computer  program  to  derive 
lower  order  systems. 

A third  order  model  provided  by  the  seeker  vendor  (Figure  12)  had  a 
slower  2.71°/sec  tracking  rate.  However  after  the  seeker  platform 
reached  the  3°  position  at  1.1  seconds,  Its  performance  was  very  close 
to  that  of  the  thirteenth  order  system  (Figure  10). 

The  state  equations  for  this  approximation  are 

FX (3)  - ,4(U700(SAMP  - X (2 ) ) - X (3) ) 

FX(2)  » X(3)  + .026316  FX (3) 

FX (1 ) - X(2) 

It  should  be  noted  that  the  state  equations  and  state  variables  for  the 
vendor  system  are  already  given  In  units  of  degrees  rather  than  In 
radian  measure  which  was  used  for  the  other  two  models. 

The  Introduction  of  both  constant  and  time-varying  drift  rates  up 
to  0.l°/sec  peak  value  had  little  effect  on  the  models'  responses. 
Greater  drift  rates,  however,  caused  failure  of  the  SAMP  forcing  func- 
tion for  the  13^  order  model  to  switch  signs  at  the  end  of  every  50  ms 
interval  (multiple  pulsing)  as  well  as  a reduction  In  the  tracking 
speeds  of  that  system.  The  vendor  model  was  not  adversely  affected 
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even  with  drift  rates  as  large  as  5°/sec.  In  actual  system  performance 
the  seeker  would  be  unable  to  cope  with  such  large  drifts,  which  Indi- 
cates the  failure  of  the  vendor  model  to  properly  simulate  internal 
friction  effects. 

Drift  was  entered  Into  the  13**1  order  system  by  the  state  equation 
in5 

FX (12)  = i ~ (-.00858  X ( 1 2 ) - X ( 1 1 ) + X(2)  + DRIFT) 

3 • / / 

and  for  the  vendor  model 

FX(2)  = X (3)  + .026316  FX(3)  + DRIFT 


The  In-house  third  order  approximation  was  not  considered  for  drifting 
since  its  response  had  so  closely  approached  that  of  the  thirteenth 
order  system. 

8.  CONCLUSIONS 


A description  of  four  computer  programs  for  the  simplification  and 
analysis  of  a discrete  data  tracking  system  has  been  presented.  These  pro 
grams  are  demonstrated  via  the  solution  of  the  step  response  of  a laser 
semi-active  seeker.  The  programs  are  general  and  should  require  minor 
modification  for  study  of  other  seeker  systems. 


FIGURE  1.  DETAILED  SERVO  BLOCK  DIAGRAM 


STATE  VARIABLE  SET-UP  OF  THE  LOOP  SERVO  SYSTEM 


oooonoooooo 


SOLUTION  Of  A DIFFERENTIAL  EQUATION  USING  THE  RUNGE-KUTTA  METHOO 

1PPS  = PULSE  REPETITION  RAIE  (PPS) 

TSAM  = SAMPLING  PERIOD  (SECONOS) 

N = ORDER  OF  DIFFERENTIAL  EQUATION 
ALOS  = ANGULAR  LINE  Of  SIGHT  ERROR  (DEGREES) 

Tn  = INITIAL  HME  (SEC) 

H = RUNGE-KUTTA  TIME  INCREMENT  (SEC) 

XN(1)...XN(N)  = STATE  VARIABLE  INITIAL  CONDITIONS  (RADIANS) 
fX (] ) ...fX(N)  = STATE  EQUATIONS  (RADIANS) 

DIMENSION  XN(20) >X(2Q) *Q(20»<?0) »F X (20) , I COUNT (20) ,XN1 (20) , Y (SOI ,5) 

68  fORMAT  <22X,6(3HXN(,  12 ♦ 1 H)  ,9X>  ,AHSAMP) 

69  fORMAT(16X,lP7EI5.6) 

?0  f OPM A T ( ] H 1 ) 

81  F ORMAT (£13) 

83  fORMAT (5E15«T) 

86  f ORMAT ( lHj 1PE1A.S* 1P7E15.6) 

Q7  fORMAT ( ///7X . ART IME • 1 OX • 7 ( 3HXN ( • 1 1 • 1H ) • I OX ) ) 

REALMS. 81  1 IPPS 
READ (S«81 IN 
READ (5«83)  ALOS 

READ (S *83)  TN,H, (XN(NN) ,NN=) .N) 

WR I TE (7.70) 

WRITE  (7,81 ) IPPS 
WRITEI7.81)  N 
WRITE (7.86*  ALOS 

C CONVERT  XN  (RADIANS)  TO  XNl  (DEGREES) 

DO  55  IQ  = I .N 
S5  XNl (IQ!  = XN(IQ)  *57,2956 

WRITE (7,86) TN.H, (XNl (NN) ,NN=l.N) 

NCOUN=0 

PPS^fLOATdPPS) 

TSAM= 1 . / ( PPS°H) 

WHITE <7.861  TSAM 

C DEFINE  THE  INPUT  f (JHC I NO  FUNCTION 

IF  < ALOS  - S7,?958*XNU  1 > 002.802.803 

802  SAMP  = -e.o 
GO  TO  BOA 

803  SAMP  =2.0 

C WRITE  THE  HEADING  AND  INITIAL  CONDITIONS 

DO  88  JX=J.N 

86  ICOUNT(JXI=lX 

WHITE' (7, 87)  (I  COUNT  (IX)  » I X = 1 .7) 

WRITE (7, 66) ( I COUNT (tX),IX=B>)3) 

WRI TE( 7,86) TN. (XNl (NN) ,NN=1 .7) 

WRITE (7,69) (XNl (Nil) ,N11=8*13) .SAMP 

C CLEAR  ThF  STATE  EQUATIONS 

6QA  DO  5 M= I iN 
S FX (N) =0.0 

C • SOLVE  THE  DIFFERENTIAL  EQUATION 

C THIS  LOADS  THE  INITIAL  CONDIIIONS  IwTO  TrtE  Y OAU  ARRAY 
Y(1  . 11=10. <»XNl  (1) 

Y ( ] • ? I =SAHP 

y<I,3)-X(3)+.5*sanj> 


FIGURE  3.  LISTING  OF  THE  RUNGE-KUTTA  SOLUTION  OF  THE 
THIRTEENTH  ORDER  MISSILE  SEEKER  SYSTEH 


r.  RUN  PROGRAM  UNTIL  TIME  TN  = H«KK«KW  SECONDS 
DO  201  KK=1 ,400 

C PROGRAM  PRINTS  EVERY  TIME  -TN  = H*KW  SECONDS 
DO  200  KW=U100 
1 L=1 
T=TN 

C SET  STATE  VARIA.L8ES  EQUAL  TO  THEIR  ITERATION  VALUES 
DO  777  K=I.N 
777  X(K)=XN(K) 

GO  TO  101 
10  00  151  K = 1 t N 
1 s 1 Q(K,L)=H«FX(K) 

T=TN*H/2. 

00  252  K=1  ,N 

252  *(k)=XN(K)*Q<K.L>/2. 

L = 2 

GO  TO  101 
20  00  ?51  K = 1 .N 
251  Q(K,l)=H*FX(K) 

T=TN*H/2. 

DO  352  K=1  ,N 

352  X(K)=XN(K) *Q(K«L)/2. 

L = 3 

GO  TO  101 
30  00  351  K=l  .N 
351  Q(K»L1=H»EX(K) 

T=TN*H 

DO  A52  K=1  ,/M 
A52  X (K)=XN(K) *Q(K,L) 
l=A 

GO  TO  101 
40  DO  A51  K=1,N 
AS  1 Q(K,LUH<*FX(K) 

GO  TO  7 

C STATE  EQUATIONS  FOLLOW 

101  FX C I 3 > - < -A  1 1 3J •SAMP)/. 006 

PX{l2)=i,E5/3,77<M-8.5ae-3«X(12)-xUl)*X(?> ) 
fx<m~-x<i2i 

TX(I0)  = 1 .E8/3»8627«(-7.7SS£-a<»X(10)-7US»>  ♦,9»1.2‘>7.9d«*  (1 1 ) ) 
FX(9I=X<10) 

Ex  «fll=I.EJ/]  .O5»(-x<0>  ♦A.03E-Qt>Fxn0)  ♦ I .?05E-6«FX  (9)  *<(9> ) 

FX (71=1 .e7/2A.502»l-5.A27E-3«X (7>-x (61-1 .6sx<s) *,aox ( 131 1 

f X (6) ~X(7> 

Ex  (S'  =J.f3/S»92«>(-X(S>  *7.7E-5,>EX(7>  *0,S2t‘-1eEx  (6)  ♦X(61) 

EX (A)  = (-X (A  I *X<51 *1.78SE“2«EX(5n /2,b 

Exrj)=2.SE3<M-X  (31  ♦6.32t'A»(-.0S5'*X(2)  *1 .52o9A.0<>(x  (A)  ♦2.3jE-2*FX 
I (Am  ) 

XX(2)5X(3>/1 »06t A 
EX  (1  lax  (23 

GO  TO  ( 10,20, JOtAOl »L 
7 TN=TN»M 
00  8 K=1*N 

6 XNCK)BXN(K»*(Q(k.n *2.«Q(K ,2) ♦2.#Q<K>31 *Q(K.A1 1/0. 


FIGURE  3.  LISTING  OF  THE  RUNGE-KUTTA  SOLUTION  OF  THE 

THIRTEENTH  OROER  MISSILE  SEEKER  SYSTEH  (Continued) 


SAMP=0. 

210  If  (NCOUN  .LT.  500)  GO  TO  200 

1F(  ALOS  - 57.2956*XN<1))  902,902*903 

902  SAMP  = -2. 0 
NCOUN=  0 

GO  TO  200 

903  SAMP  = 2.0 
NCOUN  = 0 

200  CONTINUE 

C CONVERT  XN  (RADIANS)  TO  XN1  (DEGREES) 

DO  503  IQ  = 1,N 
503  XNl(lQ)  = XN(IQ)  *57.2958 
c print  the  final  solution  of  the  state  variables 

WRITE  (7,86)  TN,  (XN1  (Nfj)  .NN=1 .7) 

WRITE <7, 69 1 (XNl (Nil ) »Nll=S,ll) »SAMP 
C This  SCALtS  THE  INPUT  OATA  ARRAY  FOR  THE  PLOT  ROUTINE 
Y (KK ♦ 1 > ) ).-10.<>XNl  ( 1 ) 
v (KK*1 j2) =SAMP 
Y(KK+1 ,3)“GUI0 

201  CONTINUE 
K-2 
NF-401 
NS-50 

WRITE(7,70) 

CALL  PLOT(V,H,NF,NS) 

STOP 

END 


FIGURE  3.  LISTING  OF  THE  RUNGE-KUTTA  SOLUTION  OF  THE 

THIRTEENTH  ORDER  HISSILE  SEEKER  SYSTEM  (Continued) 
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SUBROUTINE  PLOT  (Y.M.NF.NS) 

C V = ARRAY  OF  DATA  TO  BE  PLOTTED*  WHICH  IS  READ  IN  THE  FORM 

C Y ( 1 * 1 ) ...  Y ( NF ♦ 1 ) » Y < 1 %2>  ...  Y (NF  * 2)  ...  YU. Ml  ...  Y(NF*M) 

C M = NUMBER  Of  PLOTS 

C NF  r NUMBER  OF  PLOTTED  POINTS  ALONG  ABSCISSA.  BEGINNING  WITH  THE 
INITIAL  CONDITIONS  AT  TIME  T = 0 

NS  = MAXIMUM  UPPER  LIMIT  OF  ORDINATE  (ROUTINE  PLOTS  FROM  NS-100  TO  NS) 
DIMENSION  Y(S01  ,51  .LINE  UOl)  *LU]  ) *JL(5) 

DATA  JL  ( 1 ) * JL  (2)  * JL(  3)  ,JL  (9.)  * JL  (5)  /1HA*  1HB>  IHC,  1H0.1HF/,  JN,  JP,  J1  , 

I JBLANK. JZ/IH-. !H*,lHlilH  ,IH$/ 

00  99  1=1*101 
99  LINE ( 1 ) =UBLANK 
NN=0 
Nel 

C LABEL  THE  ORDINATE  AXIS 

DO  10]  1=1,1 1 

101  Lm  = (10»l-U0*NS) /}0. 

IF  (M-2)  103.109,109 
103  WRITE (7, \ OS)  (L(D  *1  = 1 .1 1) 

JOS  FORMAT  <2X, 11(19, 6X1, IX, 6HY(N, 1)1 

go  to  ns 

109  WRITE (7*106)  (L(I). 1=1*11) 

106  FORMAT  (2X, 11(19, 6X1, IX, 6HY(N,]>»7X»6HY'N, 2)1 
GO  TO  US 

110  IF  (NN-10)  125*115*115 
115  NO=0 

NN=  0 

00  120  1=1,10 
ND=ND* 1 
L 1 ME (ND) = JP 
00  120  0=1,9 
ND=ND‘ » 

120  LINE(NDl=JN 
LINE ( 101 >=JP 
GO  to  135 

1 25  DO  UTO  1=1*101,10 

130  LINE ( I ) = Jt 

C CHANGE  NUMERICAL  DATA  to  LETTERS 
135  DU  160  1=1 ,M 
XN$«NS 

JA  = v (N, I) * I 0 1 .99999-XNS 
IF  (JA-1011  190.155*193 
190  IF  (JA)  150*150*153 
195  LINE (101 l=JZ 
GO  TO  160 
ISO  l INC  ( ! ) =,12. 

GU  TO  160 
1SS  LINE  (JA)=JLU) 

160  CONTINUE 

C IMIS  RESCALES  Y (N,  1)  ,Y (id,?)  TO  ORIGINAL  VALUES 
Y(N,1  > = Y(N*U/10, 

NC= (N-l ) /l 0 
If  (M-2)  163*173*173 
163  IF (NN)  165.16S.17S 


FIGURE  i*.  LISTING  OF  THE  PLOTTER  PACKAGE 


165  WRITE  C7.166)  NC  • L 1 NE.  . Y (N , 1 ) 

166  FORMAT  (1A.1X.101A1.1P_E13.51-. 

go  ro  i as 

1 ?0  WRITE  (7.  171  ) LlNE.Y(N«l) 

171  FORWAT(5X,101At  .1PE1'3.51 
GO  TO  185 

173  (F  (NNl  175.175.180 

175  WRITE  (7,176)  NC  *1. 1 NE*  Y (N»  1 1 . Y (N.2 ) 

176  F ORNA  T ( 1 A . 1 X . I 0 1 A 1 . 1 P^E 1 3 • 5 1 
GO  TO  185 

100  WRITE  (7,181)  LlNE.YlN.il  .YIN.2) 

101  FORMAT(5X.10\A1  .1P2EU.51 
*85  00  190  1=1.101 

190  LINE ( 1 1 a JBLANK 
NN=NN»  t 

195  N=N» I 

IF  (N-NFl  l 10. 110.200 
200  RETURN 
END 


FIGURE  *».  LISTING  OF  THE  PLOTTER  PACKAGE  (Continued) 
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THIS  PROGRAM  SOLVES  COMPOUND  SERVO  SYSTEM  POLYNOMIALS  OF  THE  FORM 
(A!  ♦ A2  ♦ ♦ AM) (B) 


(Cl  + C2  ♦ ...  ♦ CM) <D> 

WHERE  Alt  A2»  ...  AM»  B,  Cl.  ...  CN,  0 ARE  PRODUCTS  OF  POLYNOMIALS 
AND  ANY  CONSTANT  FACTORS. 

THE  PROGRAM  FIRST  SOLVES  THE  "A"  POLYNOMIALS.  ADDS  THEM,  THEN 
MULTIPLIES  THEIR  SUM  TIMES  THE  "B'‘  POLYNOMIAL  PRODUCT. 

THE  PROCESS  IS  REPEATED  FOR  THE  DENOMINATOR 

NA  = NUMBER  OF  "A"  OR  ''CM  POLYNOMIAL  PRODUCTS  TO  BE  ADDED 
NP  = NUMBER  of  polynomials  to  be  multiplied 
KK  = NUMBER  Of  CONSTANT  FACTORS 
N1  = ORDER  OF  FIRST  POLYNOMIAL. 

N2  = ORDER  OF  SECOND  ANO  SUCCEEDING  POLYNOMIALS 

DATA  IS  READ  IN  AS  CONSTANT  COEFFICIENTS  Mi)»K<2>*  ETC.  WHICH  ARE. 
FACTORED  OUTSIDE  OF  THE  POLYNOMIALS. 

THEN  AS  (NI+l)  COEFFICIENTS  OF  THE  FIRST  POLYNOMIAL  BEGINNING  WITH 
THE  HIGHEST  ORDER.  "S'*  COEFFICIENT, 

COEFFICIENTS  OF  SECOND  OR  SUCCEEDING  POLYNOMIALS  BEGINNING  WI1H 
THE  HIGHEST  ORDER  "S"  COEFFICIENT. 

FINAL  POLYNOMIAL  IS  PRINTED  LISTING  DEGREE  OF  S AND  ITS  COEFFICIENT 

DIMENSION  A (20, 20) tPOLYl (20) »POL Y2 (20 ) »P  f 20  * » I T(20) 

REAL  KP,K(20) 

1 FORMAT (213) 

2 FORMAT! 1P8E10. 3> 

3 F ORMAT (5H  S*».I3i2X.lPEI6.6) 

4 FORMAT (19H  POLYNOMIAL  NUMBER, 13, 8H  FOLLOWS) 

5 FORMAT (IGH  PRODUCT  NUMBER, 13, 8H  FOLLOWS) 

6 F ORMAT (18H  CONSTANTS  FOLLOW) 

7 FORMAT (38H  POLYNOMIAL  NUM8ER  1 TIMES  CONSTANTS) 

8 FORMAT <//36H  FINAL  POLYNOMIAL  NUMERATOR  FOLLOWS) 

9 FORMAT <//38H  FINAL  POLYNOMIAL  DENOMINATOR  FOLLOWS) 

1 0 FORMAT  UHl  ) 

DO  300  JJJ-L,2 
READ (5, 1 ) NA 
DO  20  1=1 .20 
DO  20  J= I .20 
20  PH.JHO. 

NB  = 0 

IF(NA.EQ-O)  NA=1 
30  DO  100  JJ=1.NA 
READ (b, 1 ) NP.KK 
I POLYs I 
•'■P=i.O 

IF (KK.EQ . 0 ) GO  TO  SO 
READ  (5,2)  (KU),I-1,KK) 

Write (6,6) 

DO  AO  1=1 «KK 
KP=KP'»K(I> 

AO  WRITE  ( 6 i 2 i KM) 

50  IF(NB.EQ.l)  GO  TO  55 


FIGURE  5.  THE  POLYNOMIAL  SIMPLIFICATION  PROGRAM 
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READ (5. 1 ) Ml 
NPA1=N1 *1 

READ(5.2>  (POLYl(I),I=l,NPAn 
55  WRI TE (6»4)  1P0LY 
00  60  I = 1 «NPA  1 
ID=NPA!-I 

60  WRITE (6.3)  TO.POLYl (I) 

WRITE (6.7) 

70  DO  80  I =1 iNPAl 
I0=NPA1-I 

POl.Yl  ( I ) =KP*P0LY1  (I) 

80  WRITE'  (6.3)  IOjPOLYl  < I ) 

ISUB=NP-1 
DO  100  J=  1 • I SUB 
READ(5,1)  N2. 

NPA2=N2* 1 

READ (5.2)  f POL Y3 ( I ) *I  = 1»NPA2) 

1POL Y= I POLY  + 1 
WRITE (6j«> I POLY 
DO  90  1=1 « NP  A2 
ID=NPA2-I 

90  WRITE  ( 6> * 3 ) ID.P0LY2(  I ) 

CALL  POLSET  (NPAl*NPA2,A«POLYI*P0LYa) 

IT(JJ)=NPA1 

iprod=ipoly-i 

WRITE ( 6 » S ) IPROD 

DO  100  1=1 »NPA1 

ID=NPAl-l 

P:iD*l  .JJ)=P0LY1 (] > 

100  WRITE (6.3)  1 0 . POL  7 ) f 1 ) 
ir(NA.EO.l)  GO  TO  i-'OO 
CALL  POLAOD  (NA.P.P.POLYI.NPAI) 

NA=I 

NB=l 

GO  TO  30 

200  ir (JJJ.EQ.2)  GO  TO  210 
WRITE (6.8) 

GO  TO  220 
210  WRITE (6.9) 

220  DO  230  I=1,MPA1 
I D=NPA 1 " I 

230  WR|TF<6.3>  1D«P0LYKI) 

WRITE (6* 10) 

300  CONTINUE 
STOP 
END 

‘SUBROUTINE  POLSLT  (N1 .N2.A.POL1 .P0L2) 

DIMENSION  A<20»20) *P0LI (20 ) .P0L2 ( 20 ) 

NN=NI*N?-1 
00  11  1 = I «N2 
00  11  J=).NN 
11  A t Ji I ) =0*  0 
OO  21  I=l»Nl 
J~  1 

00  21  K*WN2 
A < J.K) =P0L1 ( l) 

21  J=J»l 

FIGURE  5.  THE  POLYNOMIAL  SIMPLIFICATION  PROGRAM  (Continued) 


00  31  1=1 *NM 

31  POL1(I)=0.0 
00  M 1=1 ,NN 
00  4J  K=1 iN2 

*1  P0L1  (I)  =POll  m+A(I.K)«P0L2(K) 

ni=nn 

return 

End 

SUBROUTINE  POLAOO  (NA,P,IT,ADD,  ICOUNT) 
DIMENSION  P ( 20 , 20 > , IT(20) , ADD (20) 
LCOUNT  = IT ( 1 ) 

DO  12  I=2»NA 

IF ( ICOUNT.GE. IT ( 1 ) > 60  TO  12 
1 COUNT= I T { I ) 

12  CONTINUE 
DO  22  1=1 .20 
22  ADD ( I ) =0. 

DO  32  I = 1 1 1 COUNT 
10= ICOUNT- 1 ♦ 1 
DO  32  J=1,NA 

32  ADDUD)=ADD<IO>»P.U*J) 

RETURN 

END 


FIGURE  5.  THE  POLYNOMIAL  SIMPLIFICATION  PROGRAM  (Continued) 
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FIGURE  6.  DATA  FOR  THE  POLYNOMIAL  SAMPLE  PROGRAM 
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NUMERATOR 
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FIGURE  7.  DATA  FOR  SIMPLIFICATION  OF  THE  TWELFTH  ORDER  POLYNOMIAL 
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OF  THE  TWELFTH  ORDER  POLYNOMIAL 
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CONTINUED  FRACTION  EXPANSION  AND  REDUCTION  IN  ORDER  OF  AN  N-TH  ORDER 
TRANSFER  FUNCTION 

NN  = TOTAL  NUMBER  of  TERMS  IN  THE  NUMERATOR 

All) ...  A fNN)  = THE  COEFFICIENTS  0 f THE  NUMERATOR  BEGINNING  WITH 
THE  CONSTANT  TERM  and  increasing  in  order 
ND  = TOTAL  NUMBER  OF  TERMS  IN  THE  DENOMINATOR 

BID... BIND)  = THE  COEFFICIENTS  OF  THE  DENOMINATOR  BEGINNING  WITH 
THE  CONSTANT  TERM  AND  INCREASING  IN  ORDER 

FINAL  REDUCED  ORDER  TRANSFER  FUNCTIONS  ARE  PRINTED  AS 
AO  + AleS  * A2*S**2  ♦ A3*$°*3  ♦...  ♦ AM^S^M 


C B<i  + B1WS  ♦ ♦...  ♦ BN*S*»N 

C 

DIMENSION  AI200) .B (200) «HI200) 

REAOI5.SOO)  NN 

500  FORMAT  1 15) 

ND=NN* 1 

HEAD  < 5 » 5 0 1 ) (AID  1 1=1  *NN> 

READ  <5*501 ) IB(I)»I-1.ND) 

501  FORMAT (AE20. 6) 

L = 0 

2 L=L ♦ l 

HIL ) -B ( 1 ) /A  1 1 > 

WRITE  1 6 .604)  L.H(L) 

60A  FORMAT(/2X,2HH(,I2,3H)  =»IPE16.6> 

A (NO) =0. 

ND=ND-1 
DO  10  1=1. ND 
B ( I ) =A ( I ) 

11*1*1 

10  A(T>-Q(I1>-A(I1>«H<L) 

1 L=L ♦ I 

H(L)=B(I  l/AU) 

WRITE (6*604)  L.H(L) 

IF(ND.EQ.l)  GO  TO  3 
N0*N0-1 
DO  20  1=1. MO 
B < I ’ =A< I > 

11=1+1 

20  A(I)=B(U)-A(I1)»H(L> 

ND=NO* 1 
B(ND)=A(ND) 

GO  TO  2 

3 MM=2*NN 

CALL  MRMUMM.H) 

STOP 

END 


FIGURE  8.  LISTING  OF  THE  TRANSFER  FUNCTION  ORDER  REDUCTION  PROGRAM 
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SUBROUTINE  MRML (N»H) 

01 MENS  I ON  HR (30,30) ,TR< 30,30) »T (30,30) ,0(30) ,TL(30«30) *C (30 ) »H(3Q) 
DO  1 I = !,N 
DO  1 J = 1»N 
1 HR ( I , J)  = 0. 

DO  10  J=l,N 
DO  20  L = 1 • J 
HR (L ,L ) = H ( J*1 -L) 

IF(L.EQ.N)  GO  TO  22 
IE ( J+l-L.EQ* 1 ) GO  TO  21 

20  HR (L ,L ♦ 1 ) = 1. 

21  J1=J+1 

DO  30  I=J1,N 
30  HR (I, I)  = 1. 

22  IF(J.GT.l)  GO  TO  23 
00  40  II  = 1>N 

DO  40  J1  = 1 ,N 

T( I ] * J! ) = HR(Il.Jl) 

40  HR ( I ] , J1 ) = 0, 

GO  TO  10 

23  CALL  MALTP  (N,N,HR,N,T • TR) 

DO  50  II  = 1 ,N 

DO  50  Jl  = 1,N 
T (11, Jl)  = TR(Il.Jl) 

50  HR ( 1 1 , Jl ) = 0. 

10  CONTINUE 

DO  110  J=1,N 
00  120  L= 1 , J 
IF(J.EQ.l)  GO  TO  121 
HR(L.L)  = H< J*l-L> 

IF( J*1-L.EQ.2)G0  TO  121 

120  HR (L*L *1 ) = 1. 

121  DO  130  I=J,N 
130  HR (I, I)  = 1. 

IFU.GT.l)  GO  TO  123 
DO  140  II  = 1 ,N 
DO  140  Jl  = l ,N 
T(Il.Jl)  = HR ( 1 1 i Jl ) 

140  HR(Il.Jl)  = 0. 

GO  TO  110 

123  CALL  MALTP  (N,N,HR,N,T,TL) 

DO  150  II  = l,N 
DO  150  Jl  = 1,N 
T(Il,jn  = TL(Il.Jl) 

ISO  HR ( 1 1 * Jl ) = 0, 

110  CONTINUE 

00  200  [=1,N,2 
L = 0 

DO  210  J=I,N 
L = L»1 

0(1)  = TR(i.J) 

210  C (L ) = TL ( I • J) 

NT  = (N-IM>'2*1 
NT1  = NT-l 
CALL  lSE(hJTl»0,C) 

200  CONTINUE 
RETURN 
END 


FIGURE  8.  LISTING  OF  THE  TRANSFER  FUNCTION  ORDER  REDUCTION  PROGRAM 
(Continued) 
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SUBROUTINE  MALTP <N,M,A»L.B,C> 

DIMENSION  A (30, 30) , B(30,30) »C<30»30) 

DO  10  1=1, N 
DO  10  J=1,L 
S=0.0 

DO  10  K=  1 ,M 
S=S*A(I.K)*B(K.J) 

10  C(I.J)=S 
RETURN 
END 

SUBROUTINE  ISE<hi,D.C) 

DIMENSION  A(30) .B< 30) *C (30) .0(30) .6(60) .DD( 30*30 ) . 
1 ON ( 30. 30) ,T(3*30> .F< 30,30) 

600  FORMAT (//////) 

601  FORMAT (1HQ.1P8E15.6) 

602  FORMAT ( 1H0, 122H  

1 - 

1 ) 

N1 =N+1 

WRI TE (6 .600) 

WRITf (6.601)  (C ( I ) . 1= 1 *N) 

20  WRI Tf (6.602) 

WRI TE (6*601)  (0(1). 1=1, Ml) 

NX=N 

DO  50  1=1,2 
DO  51  J=1.N 
NJ=NX-2*(J-1 ) 

IHNJ.IE.O)  T ( I , J)  = 0. 

51  IF(NJ.GT.O)  T(I,J)=0<NJ) 

50  NX=N+ 1 

DO  60  1=1, N 
DO  60  J=1,N 
DOIT. J) =0 , 

60  DN(X«J)=0. 

L=0 

LL=0 

LJ=i 

DO  70  1=1, N 
L=L*1 

IF (L .EQ. 3)  LJ=LJ» l 
IFtU.EQ.3)  1 = 1 
DO  71  J=L J,N 
LL=LL*1 

DD(I,J)=T(L.LU 
71  ONH  . ))=T(L,LL) 

70  LL=0 
N?=2«N 

DO  80  I=1*N2 

80  0(D=0. 

L=0 

u=o 

DO  81  1=1, N 
DO  82  J=1,N 
L=UL* J 

e?  G(L)=C( >«C< I) ♦G<L) 

81  LL=LL*1 


FIGURE  8.  LISTING  OF  THE  TRANSFER  FUNCTION  ORDER  REDUCTION  PROGRAM 
(Continued) 


31 


t 


00  90  J=] ,N 
JJ=2«(N-J) *1 
90  DN(1«J)=G(JJ) 

CALL  iNVERfDN.N.F.O.DETN) 

CALL  lNV£R(D0,N.F,0.DETD) 

XI=(-1.>**(N-1) /<2.*D(N1)>*D£TN/DETD 

RETURN 

END 

SUBROUTINE  INVER  (A.N.B.M.DET) 

OIMENSION  A(30»30)  ,B! 3000 ) .IPVOT  <30 »•  INDEX  (30.3)  .PIVOT (30) 
EQUIVALENCE  (IROV. JROW! , (ICOL.UCOL) 

57  DET-1.0 
DO  17  J=1,N 
17  IPVOT ( J) =0 
DO  1 35  1=1. N 
T=0.0 

DO  9 J=1 .N 

IF ( IPVOT ( J) -1 ) 13.9.13 
13  DO  23  K=1,N 

IF  ( IPVOT  <K ) - I ) 93.23.81 
93  IF(A8S(T)-a&S(A(U»K)) ) 83*23.23 
83  IR0W=J 
ICOL=K 
T=A ( J.K ) 

23  CONTINUE 
9 CONTINUE 

IPVOT ( I COL)  = IPVOT ( I COL ) * 1 
IF (IROW-1COL)  73,109*73 
73  0ET=-DET 
DO  12  L=1*N 
T = A(IROW,U 
A(IR0w,L)=A(IC0L.L) 
la  A( ICOL.L) =T 

IF (Ml  109*109.33 
33  DO  2 L=l.M 
T=B(IROW,L) 

8 ( I ROW.L) =8 ( 1COL  »L) 

2 8(IC0L,U=T 
109  INDEX (l.l)=IROW 
INDEX (1*21 = ICOL 
PIVOT ( I ) =A ( icou • ICOL) 

OETsDEToPIVOT(I) 

A ( ICOL  * ICOL)  = 1, 

DO  205  L=1.N 

205  A(I COL.L) =A ( ICOL .L) /PIVOT ( I) 

IF (Ml  397.397.66 
66  DO  52  L=1*M 

52  8(1C0L.L1=BUC0L*U/PIVQT(I) 

397  DO  139  LI  = 1 >N 

IF M.I-ICOL)  21.139.21 
21  T=A(LI.ICOL) 

A (LI  * ICOL) =0. 

DO  89  L=1.N 

89  A(LI*L)=A(LI*U-AUC0L.U«T 
IF(N»  139.139*18 


FIGURE  8.  LISTING  OF  THE  TRANSFER  FUNCTION  ORDER  REDUCTION  PROGRAM 
(Continued) 
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18  00  68  L=1.M 

68  B(LI,L)=B(LI.U-BUCOL.L)*T 

134  CONTINUE 

135  CONTINUE 
222  00  3 I = 1 «N 

l=N-I*l 

IF ( INDEX <L«1> -INDEX <L. 2) ) 19*3*19 

19  JROW  = INDEX (1,1) 

JCOL  = INDEX (L  »2) 

00  549  K = 1,N 


T = 


i JROW) 


A(K.JRUW)  = AIK*  JCOL) 


A(K.  JCOU 
549  CONTINUE 
j CONTINUE 
81  RETURN 
END 


= T 


FIGURE  8.  LISTING  OF  THE  TRANSFER  FUNCTION  ORDER  REDUCTION  PROGRAM 
(Continued) 
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FIGURE  9.  DATA  FOR  THE  ORDER  REDUCTION  PROGRAM 


GURE  11.  BLOCK  DIAGRAM  AMD  STATE  VARIABLE  MODEL  OF  THE 
IN-HOUSE  THIRD  ORDER  RATE  LOOP  APPROXIMATION 


FIGURE  12-  BLOCK  DIAGRAM  AND  STATE  EQUATION  MODEL  OF  THE 
VENDOR  THIRD  ORDER  SYSTEM  APPROXIMATION 


